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Abstract 


NASCAP-2K can be used to calculate plasmadynamic effects as well as steady-state 
charging and current collection, In this paper we consider electron dynamics in the 
sheath of a VLF antenna. We estimate the sheath size, and show 1-D calculations for 
both sine wave and square wave excitation. The results show strong electrostatic plasma 
oscillations at the sheath edge. Then we use NASCAP-2K to approach the problem from 
2 directions. Run as explicit Particle-in-Cell, NASCAP-2K can duplicate the square 
wave results through the first maximum in the plasma oscillation, obtaining excellent 
agreement with the 1-1) results. Due to the computational demands of PIC, this approach 
can only simulate the early transients using a common PC. To learn more about the 
floating potential of the antenna, NASCAP-2K can be run as hybrid PIC to capture the 
ion dynamics over the full wave cycle. 

Introduction 

There is current interest in generating VLF waves in space for the purpose of controlled 
experiments with the waves thought to control the trapped electron population (Inan et 
al, 2003). A transmitting antenna for this purpose would be several inches in diameter, 
tens of meters long and have bias amplitudes of hundreds of volts, perhaps kilo-volts. It 
would interact with a large volume of the surrounding plasma, and be a major driver for 
spacecraft potential. Thus, it is of major interest to be able to simulate dynamically and 
in 3-D the antenna together with its host spacecraft and the surrounding plasma, 

NASCAP-2K (NASA Air Force Spacecraft Charging Analyzer Program 2000, hereafter 
Nascap-2k) contains the Particle-in-Cell, PIC, features needed to do this type of problem. 
However, these features have not been fully exercised since the days of pulsed power 
experiment, SPEAR II (Space Power Experiments Aboard Rockets, Cohen, 1995). In 
studying the CHAWS experiment (Charging Hazards and Wake Studies, Davis, et al. 
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and space potentials, and we have developed the ability to do Hyprid-PIC, i.e., PIC ions 
with barometric electron densities. For the antenna problem we would like to follow both 
electrons and ions using PIC. It is important to note, however, that Nascap-2k currently 
only solves Poisson’s equation, rather than the full Maxwell equations, so it computes only 
curl-free, quasi-static fields. Near term development plans include electro-dynamic 
improvements to Nascap-2k. 


In this paper we describe the antenna problem and for baseline parameters, show the 
solution as calculated by a 1-D PIC code for electron dynamics in the sheath. We then 
pose a nearly identical problem for solution with Nascap-2k, demonstrate that the solu¬ 
tions agree, and show the graphics. We also use the Hyprid-PIC capability to study the 
ion dynamics for a full wave cycle. 

Statement of Problem 

The objective is to simulate the dynamic sheath around a negative thin rod. We avoid 
positive polarity because the positive half of the antenna will collect copious electrons, so 
the maximum potential it can reach is determined by numerous unknown factors, such as 
the relative size of the spacecraft and antenna. We wish to do this with realistic values of 
plasma density, ion mass, applied voltage, frequency, magnetic field, and spacecraft 
velocity. The arbitrary directions of the latter two require a three-dimensional code. A 
one-dimensional (radial) code can handle magnetic field either parallel to the antenna or 
circumferential (as would be caused by current flowing in the antenna) 

Sheath Size Estimate 


The sheath (defined as the region from which electrons are excluded) can be quite large, 
even for a fairly modest potential of about 100 volts. To calculate the sheath size, we 
specify the electric field at the antenna radius, Ro. We assume that the external space 
between the antenna and the sheath is filled with ions at ambient density, p. The electric 
field at any radius, r, between the antenna and the sheath is 


E{r) 


2 £ 0 r 


--E(Ro) 


The sheath condition is E(R s )=0, where R s is the sheath radius. We then integrate the 
electric field from R s to Ro to determine the corresponding potential. Figure l shows the 
relation between applied potential and sheath radius for a 10 cm diameter antenna. At a 
density of 10 12 m * the sheath radius at 100 V bias is about 15 cm, and grows to nearly a 
meter at a density of 10 10 m‘ 3 . The calculations to follow assume a density of 3xl0 n m' 3 , 
giving a sheath radius of about 20 cm. 
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Sheath Radius 



Figure 1. Potential vs. sheath radius for negative applied potential on a 10 cm diameter antenna. 
Curves for three differed, plasma densities are shown. 

Baseline Parameters 

Table l shows the baseline parameters for the calculation. A density of 3x10 m' is 
chosen so that the sheath is large compared with the wire but still very tractable 
computationally. In general the plasma is cold, so 0.1 eV is used in those places where 
temperature is required. The antenna frequency is set to 100 kHz, and the process is 


Table 1. Parameters for baseline calculations. 


Plasma Density 

3x10' 1 m' 3 

Electron Temperature 

0.0 or 0.1 eV 

Plasma Frequency 

3.1x10V 

Antenna Frequency 

100 kHz 

Magnetic Field 

0.0 or 0.5 gauss 

Electron Gyrofrequency 

0.0 or 8.8xl0 6 s‘ l 

Ion Species 

0 + 


followed for a half-period of 5 ps. To see the effect of magnetic field, a field of 0.5 gauss 
is chosen. The ordering of the plasma frequency, electron gyrofrequency, and applied 
frequency is co p > co t > 27# . 
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One-dimensional Calculations 

A simple one-dimensiona! finite element code was written to simulate quasistatic 
plasmadynamics about a long cylindrical antenna. The computational domain extended 
out to one meter from an antenna radius of 5 cm, and was divided into 1000 zones in 
equal increments of r 2 . Two ion macroparticles and two electron macroparticles were 
placed in each zone, with each macroparticle having equal charge. The simulation was 
run for 2000 timesteps of 2,5 ns each, making up the 5 ps half-period for the 100 kHz 
frequency. When electrons left the computational space they were replaced as thermal 
electrons at the boundary. 

Figure 2 shows the potential profile at various times in the calculation. As expected, the 
potential is rapidly screened to about the expected sheath radius as electrons are expelled 
from the sheath. At certain times a positive potential region appears. This is an effect of 
electron inertia, as the moving electrons do not stop of their own accord, but must be 
attracted back towards the sheath boundary. In Figure 3 through Figure 5 we plot (1) the 
maximum potential at times when a positive region appears; (2) the location of the 
maximum potential; and (3) the location of the sheath edge, indicated by a sharp drop in 
the charge density front nearly the ambient ion density to nearly zero. 



0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.65 0.95 


Radius (m) 

Figure 2, Potential profi e at various times during the one-dimenional calculation. 
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Figure 3. Simulation results for half sine wave and no magnetic field, showing peak positive 
potential (magenta curve, right scale), location of peak (yellow curve, left scale) and location of 
sheath edge (dark curve, left scale). 



Tlmi (a) 


Figure 4. Same as Figure 3, for half sine wave with magnetic field of 0.5 gauss. 

Figure 3 and Figure 4 show results for a half sine wave, for which the applied negative 
potential continuously rises and returns to zero. The magnitude of the potential maximum 
is 8 to 10 V, and the oscillation frequency is somewhat less than the electron plasma 
frequency. The sheath edge occurs at a radius of about 20 cm as calculated above, with 
oscillations of about two cm. The potential maximum, when it occurs, is just inside the 
sheath edge. 
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Figure 5. Same as Figure 3, for square wave with 0.5 gauss magnetic field. 

Two differences can be noted between the un-magnetized (Figure 3) and magnetized 
cases. First, in the un-magnetized case the potential goes completely non-positive 
between peaks, whereas in the magnetized case a positive peak forms inside the outer 
boundary, This occurs because the magnetic field inhibits inward diffusion of the thermal 
electrons. Second, in the un-magnetized case the sheath remains at the end of the pulse, 
even though the potential goes to zero, while in the magnetized case the sheath 
disappears. This occurs because the magnetic field traps inbound electrons in the sheath 
region, whereas in the absence of magnetic field inbound electrons collide with the 
antenna. 


Square wave calculations were done in one-dimension to compare with the Nascap-2k 
calculations below, and are shown in Figure 5. Much stronger plasma oscillations are 
seen, with the initial oscillation at 75 V. The sheath edge oscillates 10 cm on either side 
of its average position at 20 cm, and the peak potential occurs well inside the sheath. The 
0.5 gauss magnetic fie d doubles the rate of decay of the oscillations compared to the un- 
magnetized case (not shown). 

Three-dimensional PIC Calculations 

Three-dimensional calculations were done with Nascap-2k to demonstrate the feasibility 
of such calculations. Figure 6 shows the Nascap-2k antenna model embedded in a nested 
grid. The antenna consists of two square rods, each 10 cm on a side and 4 m long. The 
outer boundary of the grid is a square 1.32 m on a side. The coarse resolution is 11 cm, 
with 5.5 cm resolution near most of the antenna, and 2.75 cm resolution in a limited 
region. Initially, 8 electron macro-particles and 8 ion macro-particles were placed in each 
zone, positioned so as to represent a uniform charge distribution in the context of the non¬ 
linear interpolants. A negative 100 V square wave was applied to half of the antenna, and 
each timeslep consisted of (1) tracking the particles for 2.5 ns, (2) sharing the particle 
charge to the nodal coefficients in accordance with the nonlinear interpolants, and (3) 
recalculating the potential in preparation for the next tracking phase. 
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Figure 6. Nascap-2k antenna model, showing antenna and gridding. 

The initial state is represented in a three-dimensional view in Figure 7, and a planar view 
in Figure 8. These figures show the potential in a cut plane, and a plane of nearby 
electron macro-particles. The ion macro-particles have the same initial configuration, but 
move negligibly during the simulation time. Note that the apparent high density of 
particles in the subdivided region is balanced by correspondingly reduced particle weight. 
Figure 9 show close-ups of the simulation at later times. The left figure shows the 
electron motion at 50 ns leading to sheath radii of 18 cm. Note the effects of the square 
cross-section. Particles that started out near the flat, low-field region have moved 
considerably less than those that started out near the high-field comers. Figure 9-right 
shows the configuration at time of maximum positive potential (137 ns). At this time, 
there is a high, broad maximum in the potential, with electrons excluded from a region 
that extends well beyond the location of the potential maximum. 
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Figure 7. Nascap-2k antenna model showing potentials and particle positions after 2.5 ns. 


Figure 8, Planar view of initial potentials and particles, as shown in Figure 7. 
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Figure 9. Blowups of particles and potentials. Left: sheath radius of about 18 cm 
after 50 ns. Right: at time of maximum positive potential (137 ns). 

Figure 10 shows another view of the final configuration, with potentials in a plane 
containing the antenna Note that there is no apparent difference between the potentials 
the highly resolved region and in the less resolved region, suggesting that the highest 
level of resolution ma) not be needed. 



Figure 10. Another view of the final configuration (at 137.5 ns), showing potentials in a plane 
containing the antenna. 
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Figure 11 shows the development of the sheath radius and maximum potential calculated 
by Nascap-2k, compared with the one-dimensional sheath radius result. The maximum 
potential of 75 V, as well as the time of first maximum (140 ns) is in excellent agreement 
with the one dimensional result. The maximum sheath radius calculated by Nascap-2k is 
larger than the one-dimensional result in proportion to the effective larger size of the 10 
cm square antenna vs. the 10 cm diameter round antenna. 



Time (ni) 


Figure 11 . Nascap-2k results for sheath radius (dark curve) and maximum potential (magenta 
curve, right scale) compared with one-dimensional sheath radius results (yellow curve), 

Hybrid-PIC Calculation 

We also attempted a more ambitious antenna calculation than above. The intent was to 
show that a calculation could be performed for parameters more nearly resembling the 
frequency, waveform, and plasma density of interest for VLF experiments. The old and 
new calculations are contrasted in Table 2. In both cases, the spacecraft body and positive 
portion of the antenna were grounded, as static calculations indicated that electron 
currents would limit the body to negligible positive potentials. 

The Hybrid PIC technique is appropriate for low frequency problems where the electrons 
may be expected to behave as an equilibrium fluid, but where the ions must be treated as 
a kinetic gas using the PIC method. In this case, we used the simplest electron fluid, the 
Boltzmann equilibrium, and only the 0 + ions were tracked. The potentials were again 
solved at each time step. The calculations were carried out to 1.9 cycles of the wave. A 
snapshot of the 2 kHz calculation is shown in Figure 12. The figure shows the potentials in 
a cut plane midway through the negative end of the antenna. Although only a 
representative subset of the ions macro-particles are shown, the figure still illustrate how 
Nascap-2k can perform these calculation with a relatively low number of particles 
compared to more common approaches to PIC due to the high order finite element potential 
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solver. The Nascap-2k potential solver uses non-linear interpolating functions that can 
guarantee strictly continuous electric fields as well as potential. This dramatically reduces 
the electric field noise that would otherwise force the calculation to use orders of 
magnitude more particles and accordingly more times. These simulations took on the order 
of 1 day on a dedicated GHz PC. 


Table 2. Contrasting parameters for the old and new calculations. 


Parameter 

PIC Calculation 

Hybrid-PIC Calculation 

Frequency 

100 kHz 

2, 6. 20 kHz 

Plasma Density 

3x10" m' 3 

IxIO 10 m' 3 

Plasma Frequency (electron) 

4.92 MHz 

898 kHz 

Plasma Frequency (Oxygen) 

28.5 kHz 

5.2 kHz 

Peak applied poiential 

-100 V 

-100 V 

Waveform 

Square wave 

Sine Wave 

Sheath radius 

34 cm 

60 cm 

Timestep 

2.5 ns 

l/50f 

Time Simulated 

0.140 ns 

2/f 

Antenna Length 

4 m 

10 m 

Antenna Radius 

0.05 m 

0.1 m 

Initial No. of Electrons (Ions) 

350,000 

966,000 



Potentials 

.VBM- 

5.0 

..I 


Figure 12. Subset of ion positions and potentials from 2 kHz model at time step 20 at an antenna 
voltage of -59 Volts. 


The calculations do noi reached a fully repeating state, but the cycle to cycle variations are 
not so significant that we cannot draw some preliminary conclusions. Figure 13 shows a 
summary of the three calculations at frequencies of 2, 6, and 20 kHz. The ion plasma 
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frequency is 5 kHz. It can be observed that the current lags the voltage by an amount that 
becomes noticeably significant as the frequency exceeds the ion plasma frequency. We can 
also observe that ion current profile begins to deviate from the driving sinusoid and appears 
to be developing harmonics close to the ion plasma frequency. 
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Figure 13. Current and voltage history for the Hybrid PIC calculation. 

Conclusions 

When the driving waveform contains high frequency transients as with the square wave, 
electron dynamics are important in the large sheaths of VLF antennas under space 
conditions. Large electrostatic oscillations are to be expected, producing at times positive 
potential regions aboul a negative antenna. The amplitude of the oscillations is waveform 
dependent. Nascap-2k is shown to have the ability to perform quasi-electrostatic PIC 
simulations of sheath dynamics in fully three-dimensional geometry. This is important 
because such simulations can take into account the presence of the host spacecraft and 
arbitrary directions of magnetic field and spacecraft velocity. 

For waveforms with only low frequency content, we have assumed that the electron 
dynamics can be ignored as a first approximation and have shown that the ion sheath also 
has interesting dynamics features. These are very preliminary calculations but a tentative 
conclusion seems to be that without some means to hold the antenna at a more positive 
potential, the negative potentials will drive ion dynamics that may be rich in harmonic 
content. Ultimately, predicting the radiation of a low frequency antenna will require 
following both the electron and ion dynamics as well as a more electro-dynamic 
treatment of Maxwell’s equations. This will drive calculations times, which have been 
10s of hours on fast PCs for these simulations, upward by a factor over 100 unless other 
innovations are introduced. Some of the techniques under evaluation are orbit averaging, 
where a particle’s charge is spread over the full path covered by a large trajectory step; 
gyro-averaging for the electron orbits, and of course parallel processing. 
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